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The soft-core boson system is one of the simplest models of supersolids, which have both off- 
diagonal long-range order (Bose-Einstein condensation) and diagonal long-range order (crystalline 
order). Although this model has been studied from various points of view, studies of the stability 
of current-flowing states are lacking. Solving the Gross-Pitaevskii and Bogoliubov equations, we 
obtain excitation spectra in superfluid, supersolid, and stripe phases. On the basis of the results of 
the excitation spectra, we present a stability phase diagram that shows the region of the metastable 
superflow states for each phase. 
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Introduction. A supersolid is a quantum phase that 
has both superfluidity and solidity. After it was demon- 
strated in some seminal works 1 - - — that this intriguing 
state can be realized in quantum solids, much research 
has been done. Recently, Kim and Chani presented the 
possibility of a supersolid in solid 4 He. They reported 
experimental results of non-classical rotational inertia 
(NCRI) in solid 4 He. While there have been many ex- 
perimental and theoretical studies made££, there still re- 
mains controversy over the origin of NCRI in solid 4 He. 
Ultra-cold atomic gases have also become a field for the 
study of supersolids. Many theoretical works support the 
existence of a supersolid phase in systems with strong 
long-range interactions such as dipole-dipole^ or van der 
Waals interactions^. Recently, Bose-Einstein conden- 
sates (BEC) of 52 Cr£, 164 Dyi2, and 168 Erii, which have 
large magnetic moments, have been realized experimen- 
tally. Owing to the high controllability of experimental 
conditions, cold atomic gases may provide insights into 
the nature of supersolids. 

In the many previous studies, the properties of equilib- 
rium states of supersolids have been investigated. How- 
ever, dynamical stability of a finite superflow state can 
not be obtained by calculating equilibrium states in the 
absence of a current 1 ^. We need to investigate the stabil- 
ity of a current-carrying state that has both off-diagonal 
long-range order (ODLRO) and diagonal long-range or- 
der (DLRO). 

The stability of condensates can be investigated by 
calculating the excitation spectra. In many cases, the 
low-energy excitations determine the critical velocity of 
supcrfluids. For example, the Landau critical velocity^ 
is determined by the excitation energy, and the critical 
velocity of a condensate in a moving optical lattice can 
be calculated from the excitation spectrum 1 ^. The ex- 
citation spectrum can be not only theoretically calcu- 
lated but also experimentally observed in neutron scat- 
tering experiments for 4 Hei£ or Bragg spectroscopy of 
cold atoms 1 ^. One of the most striking properties of su- 
pcrfluids is the existence of a critical velocity, and deter- 
mining the critical velocity of a supersolid is an impor- 
tant problem. To determine the supersolid critical veloc- 



ity, we use a simple continuum model of supersolids that 
was proposed by Pomeau and Rica-i. Using the Gross- 
Pitaevskii (GP) equation with a finite-range interaction 
(soft-core interaction) , they showed that the ground state 
exhibits ODLRO and DLRO at a sufficiently high density 
of the condensate. Although this model uses a simplified 
interaction potential for the purposes of manageability, 
it is suitable for investigating supersolidity. In fact, the 
properties of bosons with finite-range interactions have 
been studied in various contexts ^ 18 ' 19 ' 21- — . 

In this paper, we present a phase diagram of 
metastable superflow states for each phase(we call this 
the "stability phase diagram" in the following) of two- 
dimensional soft-core bosons at zero temperature by solv- 
ing the GP and the Bogoliubov equations^. Although 
a similar analysis has been done for a lattice system 2 ^, 
a continuum system has not yet been studied. We find 
that three phases are stable against the superflow: a su- 
perfluid phase, a supersolid phase, and a stripe phase. 

Model. We use the two-dimensional GP equation with 
a finite range interaction 2 *^ -1 ^ 21 : 

V 2 *(r) + dr'V(r - r')|*(r')| 2 vP(r) = ^(rll) 

2m J 

where ^(r) is the condensate wave function, m is the 
atomic mass, V(r — r') = Vo0(a — \r — r'\) is the two- 
body interaction, Vq is a positive constant, a is the inter- 
action range, and 9(x) denotes the Hcaviside step func- 
tion. The chemical potential /i is determined by the to- 
tal particle number TV. The interaction strength of this 
system can be measured by a dimensionless parameter— 
g = nnoma 4 Vo / h 2 , where no is the mean-particle den- 
sity. We use j as a control parameter. We assume that 
the solution of eq. ([TJ is a plane wave or has a crystalline 
order that can be written by the Bloch wave function^: 



*(r) 



q+Gt 



(2) 



G 



where fiq/m is the velocity of the condensate, G is the 
reciprocal lattice vector, and C q +G is an expansion co- 
efficient. We calculate three crystalline structures: a tri- 
angular lattice, a square lattice, and a stripe structure. 
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We optimize the period and shape of these structures to 
minimize the energy per particle under given parameters 
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q and g. In order to obtain the excitation spectrum, we 
solve the Bogoliubov equation: 



eu(r) = Ku(r) + / dr'V(r - r') [**(r')*(r)w(r') - *(r')*(r)w(r')] , 



ev(r) = -Kv{r) - / dr'V(r - r') [#(r')**(r)v(r') - #*(r')**(r)u(r')] , 



K 



2m 



V 2 -/i + J dr'V(r-r')\V{r')\ 2 , 



(3) 
(4) 
(5) 



where e is the excitation energy and u(r) and v(r) are 
excitation wave functions. Using the Bloch theorem, we 
can expand u(r) and v(r) in terms of the reciprocal lat- 
tice vector: 



Uk,n{r) = 
Vk,n{r) 



i{k+G)-r 



G 



-iq r q 



k+G,n£ 



i(k-G)-r 



(6) 



(7) 



G 



where k is the wave number vector of the excitations, 
n is the band index, and Ak+G.n and Bk+G,n are ex- 
pansion coefficients. Solving the GP and the Bogoliubov 
equations for the four assumed states, we obtain the ex- 
citation energy and local stability of each phased. We 
restrict the number of expansion coefficients to 29, 73, 
and 89 for the stripe structure, triangular lattice, and 
square lattice around the origin of the reciprocal lattice 
space, respectively. We have checked the convergence of 
the present results by comparing the results with 27, 61, 
and 81 expansion coefficients for the stripe structure, tri- 
angular lattice, and square lattice, respectively 

Results. First, we show the ground state, which cor- 
responds to the case of q = 0. In our calculation, we 
find that a square lattice structure is always dynamically 
unstable and the stripe phase is not realized at q = 0. 
Henceforth, we use the term "supersolid(SS) phase" only 
for a triangular lattice. The definitions of these phases 
are summarized as follows: the supcrfluid(SF) phase is 
the state in which the density of the condensate is uni- 
form, the SS phase has a triangular lattice structure, and 
the stripe phase has a one-dimensional periodic struc- 
ture. A typical density profile of the SS phase is shown 
in Fig. [T] (a). Comparing the energies of the stable sta- 
tionary solutions for SF and SS phases, we find that the 
SF (SS) phase has a lower energy than SS (SF) phase 
at g < g c ~ 39.49 (g > g c ). However, the chemi- 
cal potentials of the SF and SS phases are not equal 
at g = g c . This implies that an inhomogencous phase 
( = coexistence phase) is realized as the ground state 
in the vicinity of g c . We determine the coexistence re- 
gion as 38.44 < g < 40.98; the SF phase is realized for 
< g < 38.44 and the SS phase for 40.98 < g. The 
result of the transition point is consistent with that of 
Ref.|2j&. 




FIG. 1. (Color online) (a) Density profiles of the SS phase 
for (g,qa) — (45,0) and (b) the stripe phase for (g,qa) — 
(37.5,1.05). 



Next, we consider the current-flowing states with q ^ 
0. The current is assumed to be parallel to the x di- 
rection (q = (q,0),q > 0). We do not consider coexis- 
tence phase but single phases in the following (we will 
discuss this point later). Figures [2] and [3] show the sta- 
bility phase diagram that represents the region of single 
phase of metastablc superflow states in the (g, qa) plane. 

In the SF phase, the condensate wave function and 
the chemical potential are given by \&(r) = y/noe lq ' T and 
/.t = TTnoVoa 2 + h 2 q 2 /(2m). Substituting these expres- 
sions into the Bogoliubov equation, we can obtain the 
analytical expression of the excitation spectrum in the 
SF phase: 



<k 



h 2 q ■ k 



'h 2 k 2 
2m 



n 2 k 2 

2m 



ATrn Voa 2 



J\ (ka) 
ka 



,(8) 



where J\{x) is a Bessel function. In the case of q = 0, 
this spectrum has a roton minimum when g > 15.81 and 
the roton gap vanishes at g ~ 46.30. The metastable 
region of the SF phase is bounded by the thin solid red 
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FIG. 2. (Color online) Stability phase diagram in the (g, qa) 
plane. The thin solid red , dashed blue, thick solid green, dot- 
ted blue, and dashed-dotted green lines represent the Landau 
critical velocity of the SF phase, Landau instability(LI) line of 
the SS phase, dynamical instability(DI) line of the SS phase, 
DI at the long wavelength limit line of the stripe phase, and DI 
at finite k line of the stripe phase, respectively. The SF phase 
is metastable in the region surrounded by the thin solid red 
line and the line with qa — 0. The SS phase is metastable on 
the right side of the dashed blue line. There is no metastable 
stationary solution in the region not surrounded by any lines. 
Since we can not see the metastable region of the stripe phase 
in the present plot range, the magnified figure is shown in 
Fig. I 
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FIG. 3. (Color online) Magnified view of the region near the 
stripe phase of Fig. [2] The stripe phase is metastable between 
the thin solid red line and the dotted blue and dashed-dotted 
green lines. 



lines in Figs. [2] and [3l which represent the Landau crit- 
ical velocity^. Outside of this boundary, SF phase is 
unstable against the supcrflow. 

The metastability of the SS phase can be judged from 
the excitation spectra which are obtained by numer- 
ical calculations. Figure [4] (a) shows the typical excita- 
tion spectrum of the SS phase at q = 0. There are three 
branches of gapless excitation in the long wavelength 
limit. These modes are Nambu-Goldstonc modes: one is 



a Bogoliubov mode due to spontaneous U(l) symmetry 
breaking and the others are transverse and longitudinal 
phonon modes due to spontaneous breaking of transla- 
tional symmetry . In the absence of a supercurrent, the 
excitation spectrum in the SS phase of this model and 
that of more general systems were studied in Refs. [U 
and [3(| respectively. We note that the lowest branch 
in the SS phase is the Bogoliubov mode, which causes 
instabilities in the SS phase as described later. 

The metastable region of the SS phase is shown as the 
region on the right side of the dashed blue line or thick 
solid green line in Fig. (3p2. Neither stable nor metastable 
stationary solution of the SS phase exists on the left side 
of these lines. The dashed blue and thick solid green 
line represent the Landau instability(LI) and dynamical 
instability(DI) lines of the SS phase, respectively; LI and 
DI, respectively, denote the instabilities caused by the 
excitation whose spectra (ek,n) have negative real part 
and nonzero imaginary part. Figure 4 (b) shows that the 
Bogoliubov mode that has a negative real part, which 
destabilizes the SS phase. These instabilities (LI and 
DI) occur at the long wavelength limit in the SS phase. 

The metastable stripe phase exists in the region sur- 
rounded by the thin solid red line and the dotted blue 
and dashed-dotted green lines in Fig. [3] Typical den- 
sity profile of the metastable stripe phase is shown in 
Fig.Q](b). The dotted blue and dashed-dotted green lines 
in Fig. [3] represent DI at the long- wavelength limit and 
DI at finite fc, respectively. Both of these instabilities oc- 
cur in the longitudinal phonon mode unlike in the case of 
the SS phase. The mechanism through which the stripe 
phase is generated was discussed in Refs. l3lT[33l : namely, 
a negative excitation mode occurs at finite fc, leading to 
instability and the growth of a one-dimensional periodic 
structure along the flow direction. This is the reason why 
the stripe phase exists above the LI line of the SF phase. 
The period of the stripe phase is determined by the wave 
number that yields the roton minimum, as expected from 
Ref. l32l . The excitation spectrum in the stripe phase is 
plotted in Fig. [5] In contrast to the SS phase, there are 
two gapless modes in the long wavelength limit. The low- 
est gapless mode is the longitudinal phonon mode and 
the other is the Bogoliubov mode. Since translational 
symmetry is broken in only one direction, the only two 
gapless modes exist in the long wavelength limit. 

Although we have focused on the metastable states 
in the form of Eq. ([2]), there are many branches of the 
metastable states that cannot be written by Eq. @ such 
as a coexistence phase for nonzero q. Further, stationary 
solutions including a point defect or a complex network 
of defects have been reported in Ref. 0- In reality, which 
states are realized among metastable states depends on 
the initial condition and experimental procedures, for ex- 
ample, how the velocity of a container is developed from 
zero as a function of time. In order to determine the 
final state, we need to calculate a real-time and real- 
space dynamics with a specified protocol. This is a future 
work. Our results, on the other hand, implies that we can 
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reach long-lived supcrfluid, stripe and supcrsolid phases 
only when the values of {g, qa) correspond to the regions 
shown in Figs. 2 and 3 for respective phases. Particu- 
larly, we see that realization of a stripe phase requires a 
fine tuning of the parameter (g, qa) in the present model. 
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FIG. 4. (Color online) (a) Excitation spectrum of the SS 
phase for (g,qa) = (45,0). eo = Tt 2 /ma 2 is the energy 
unit. The coordinates of T, M, and K points are given by 
(0,0), (27r/v / 3A,0), and (2tt/v / 3A, 2tt/3A), where A is a lat- 
tice constant. Similar results have been shown in Ref. l2ll . 
(b) Excitation spectrum of the SS phase along k y — for 
{g,qa) = (45,1.17). 




FIG. 5. (Color online) Excitation spectrum along k y = in 
the stripe phase for (g,qa) = (37.5, 1.05). 

Summary and discussion. In summary, we investigated 
the nature of the two-dimensional soft-core bosons at zero 
temperature by solving the GP and Bogoliubov equa- 
tions. The superfluid, supersolid, and coexistence phases 
appear as the ground states. We presented the stability 
phase diagram, which represents the region of the homo- 
geneous metastable states. The metastable superfluid, 
supcrsolid, and stripe phases are realized. 

The problem of the presence of impurities or obstacles 
in the supersolid phase is a consideration for future work. 
The authors of Ref. [l?] concluded that the superfluidity 
of a supersolid phase in the presence of an obstacle van- 
ishes. However, in our calculation, a supersolid phase still 
exhibits a supercurrent in the presence of an obstacle^. 
This discrepancy could be solved from the viewpoint of 
the stability analysis. 
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SUPPLEMENTAL MATERIAL 



I. NUMERICAL METHOD FOR THE GP EQUATION 



In this appendix, we show the numerical method for the Gross-Pitaevskii(GP) equation. A similar method was 



Dpei 

used in Ref. Ill In the following, we consider two-dimensional system. 
The time-dependent GP equation is given by 



dt v 1 



2m 



V 2 *(r,f) + / dr'V(r - r')|*(r' ,t)| 2 *(r,i). 



(9) 



(10) 



Substituting *(r,t) = e~ ilJtt / h ^(r) into Eq. ©, we obtain the time-independent GP equation 

fi 2 f 

V 2 *(r)+ / dr'V(r-r'm(r')\ 2 ^(r) = ^(r), 

2m J 

where /i is a chemical potential determined by the total particle number condition 

N = J dr\V{r)\ 2 , (11) 

where N is the total particle number. We assume that the ground state wave function can be written by the Bloch 
wave function : 



*(r) = e iq - r (j,(r) 



(12) 



where <fi(r) is a periodic function that satisfies 4>(r + tiiai + 71202) = 4>{ r ) f° r arbitrary integers n\ and 712 and 
primitive vectors a\ and a 2 . From this assumption, we can expand 'l'(r): 



(13) 



where G is a reciprocal lattice vector and C q +G is a expansion coefficient. In the following, we abbreviate C q +c to 
Cg- Substituting Eq. (fl"3|) into the GP equation and multiplying both sides of the resultant equation by / dre~ lG r , 
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we obtain 



2m 



C G + J2 SagCg+ag = »Cg, (14) 

AG^O 

S AG = V(AG) £ C G , +AG C G > , (15) 

G' 

where we introduce the Fourier transform of the two-body interaction 

V(k) = J dre~ lk r V{r), (16) 
and we used the orthonormality and the total particle number condition 

-/ dre^ G - G '> r = S G ,G', (17) 

S Ju.c. 

«o^«£|C G | 2 , (18) 

G 

where S is the area of the unit cell and / dr denotes the integral over the unit cell and no is the mean-particle 

Ju.c. _ 

density. In the case of the soft-core interaction, V(k) is given by 2 - 

V{k)=2*V«a 2 '±p±, (19) 
ka 



where J\{x) is a Besscl function. For example, S is given by 



S = X (square lattice), (20) 
\/3 

S = ~^~^ 2 (triangular lattice), (21) 



where A is a lattice constant that is determined by minimizing the total energy per particle: 

■ f dr|W(r)| 2 + -L f dr f dr'V(r - r')|*(r')| 2 |^(r)| 2 

J u.c. J u.c. J 



E h 2 



N 2mN 



h 2 „ l 



5> + G) 2 |C G | 2 + — V(Gi-G 3 )C- Gl+G2 _ G3 C G3 C G2 C Gl . (22) 



G G1.G2.G3 



Here, regarding the GP equation (TH| as the eigenvalue equation for /i, we solve the GP equation numerically. Since 
the GP equation is the non-linear equation for Cg, we need to solve it self-consistcntly. The advantage of this method 
is to apply the same diagonalization algorithm to solving the Bogoliubov equation described later. The detail of the 
procedure is as follows: 



(i 
(ii 
(m 
(iv 
(v 
(vi 
(vii 
(viii 



Choose appropriate value of A. 

Choose appropriate {Cg} as an initial condition. 

Substitute {Cg} into Sag in the left-hand side of Eq. (fT4|) and diagonalize Eq. (|14[) numerically. 
Calculate the total particle energy per particle for each eigenstate. 
Choose {Cg} for the lowest energy state. 
Iterate (iii)-(v) until convergence. 

Choose different value of A to minimize the total energy per particle. 
Iterate (ii)-(vii) until convergence. 
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II. NUMERICAL METHOD FOR THE BOGOLIUBOV EQUATION 

In this appendix, we show the numerical method for solving the Bogoliubov equation. Substituting 

*(r,t) = e-*"*/" h(r) + u(r) e -' iet / h - v*(r)e iet/h ] (23) 
into Eq. ([9]) and retaining u(r) and v(r) up to 0(u(r)) and 0(ii(r)), we obtain the Bogoliubov equation 

eu(r) = Ku{r) + [ dr'V(r - r') [**(r')*(r)u(r') - $>(r')$(r)v(r')] , (24) 



ev(r) = -Kv(r) - / dr'V(r - r') [$(r')**(r)v(r') - **(r / )1'*(r)u(r')] , (25) 



2m 



V 2 - n + J dr'V(r-r')\^(r')\ 2 , (26) 



Assuming the crystalline order, we can apply the Bloch theorem to the Bogoliubov equation. u(r) and v(r) can be 
expanded by the reciprocal lattice vector— : 



l (r)=e^ r J2^ +G e z{k+G) - r , (27) 



G 



v k , n (r) = e- lq r B k+G e^ k - G ^ r , (28) 



G 



where Tik is a quasi-momentum of the excitations and n is the band index. In the following, we abbreviate A k + G and 
Bk+G to Aq and Bg, respectively. Substituting Eqs. (|2T|) and (|28|l into the Bogoliubov equation and multiplying 

both sides of Eqs. (JUJ) and (|25) by / dre~ ig r e- tfc - , 'e- tG r and / dre t<3 r e - lfc r e tG r , respectively, we obtain the 

.V u.c. J u.c. 

Bogoliubov equation for the expansion coefficients : 

D+A G + S G,AG A G+AG-J2 W G.AG B G+AG = ek,nA G , (29) 

AG#0 AG 

-D G B G - J2 S G:AG B G+AG + Y, W G*AG A G+AG = e Kn B G , (30) 
AG#0 AG 

# G = ^(q±k + Gf - fi + noV(0)+Y / V(±k + G-G')\C G '\ 2 , (31) 

G' 

'+agCg', (32) 

G' 

^G,AG = E ^(± fc + G - G')C 2G+ AG-G< C G , . (33) 
G' 

Substituting the solution of the GP equation into Eqs. ((29)) and (|30|) and diagonalizing these equations, we obtain the 
excitation spectrum 
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